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We extend the methods recently introduced m Phys. Rev. Lett. 106 036401 (2011) to investigate 
correlations between two spin-polarons in a quasi-two-dimensional Cu02 layer. The low-energy 
wavefunctions for two doped holes introduced in a half-filled Cu02 plane with 32 copper and 64 
oxygen sites are calculated explicitly using an efficient yet accurate truncation scheme to model 
the antiferromagnet background. The energetics and wavefucntions show that the charges form 
three-spin polarons and the spin is carried by a disturbance around the three-spin polaron core. 
The low-energy band results from the competition between the kinetic energy and a local attractive 
potential which favors d^2_y2 states. Lastly, we point out features that are expected to be robust 
for larger systems. 
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I. INTRODUCTION 

Spin-i antifcrromagnctism is of great importance in 
the quantum description of many exotic materials. Set- 
ting as a reference point the exhaustively studied nearest- 
neighbor Heisenberg antiferromagnet (AFM) on the two- 
dimensional square lattice^^ major developments include 
studies of frustration due to lattice geometry and to 
long(er)-range coupling, as well as the disturbance due 
to additional fermionic charge carriers. Investigations of 
the latter problem have been strongly driven by the need 
to understand the doping-controUed evolution of cuprate 
layers. This is also the focus of this article. 

Since their discovery in 1986^ cuprates have been clas- 
sified as high-Tc superconductors. Their main challenge 
to condensed matter physics is to understand the ba- 
sic mechanism leading to high temperature superconduc- 
tivity, but equally important is the need to understand 
the many anomalous properties when these materials are 
tuned away from the superconducting phase. The class of 
hole-doped cuprates, which allow electron removal from 
the parent compound, is of particular interest because of 
the clear separation of the pseudogap regime, in addition 
to the antiferromagnetism, superconductivity, Fermi liq- 
uid and non-Fermi liquid phases which occupy different 
regions of the phase diagram. The connections between 
these different phases have not yet been fully elucidated 
despite many proposals^"— In fact, one of the few widely 
accepted ideas in this field is that to find the pairing 
mechanism will require understanding the various non- 
superconducting phases first. 

It is also widely believed that a complete description of 
the lightly hole-doped spin-i 2D antiferromagnet (AFM) 
with full quantum fluctuations could provide clues for un- 
derstanding the origin of the non-Fermi-liquid behavior 
and the superconducting ground state observed at higher 



the few holes-f-AFM problem is a crucial first step needed 
to appreciate the significance and importance of such ad- 
ditions. This problem is difficult because of the compli- 
cated nature of the 2D AFM background, whose quantum 
fiuctuations in the presence of multiple doping holes were 
never fully captured for a large Cu02 lattice. 

Such a theoretical or computational description is chal- 
lenging because of the strongly-correlated many-body na- 
ture of the system. Microscopic hole- AFM interactions 
have been studied in models with oner^"— two^^"— 
three r^^— or more^^Ti^ bands. While exact analytical so- 
lutions seem to be out of reach, numerical studies are also 
carried out with various compromises, such as the use 
of small clusters, variational approaches, and/or model- 
ing of the AFM state as a classical Neel state plus spin- 
wavesi^"— Due to these limitations, there are uncertain- 
ties regarding the minimal microscopic model; after all, 
it took decades to gain a satisfactory understanding of 
even the simplest one-band t-t'-J model.—"— While cer- 
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doping. Consideration of more exotic scenarios 
exciting developments; however, a detailed modeling of 



Loss Spectroscopy, Scanning Tunneling Microscopy, re- 
cent neutron scattering and x-ray absorption measure- 
ments cannot be described by one band models, the sig- 
nificance of omitting other bands cannot be quantified 
without access to solutions of more detailed models. 

Cuprates exhibit charge-transfer band-gap behavior 
with mobile holes located mainly on the anion 2p or- 
bitals and unpaired electrons on cation 3d orbitalsi^ 
One-band models use superexchango^ and Zhang-Rice 
singlets (ZRS)^ to reduce the {N — n)-electron prob- 
lem to one of n holes in an AFM background, in which 
both spin and charge degrees of freedom are assumed to 
be hosted in the same single band. To reach agreement 
with experiments, such models must be tweaked at least 
by adding longer-range hopping^^"— i^i^ and possibly by 
coupling to phononsiii^i^. 

We have recently shown that even for the one-hole case, 
distinguishing cation and anion sites leads to significantly 





FIG. 1. (color online), (a) Two unit cells of the Cu02 plane. 
The orbitals kept in the three-band model of Eq. ^ are 
shown, with white (shaded) for positive (negative) signs. The 
two e vectors (solid arrow) and four 5 vectors (dashed arrow) 
are also shown, (b) Sketch of a virtual process of Tswap. 



different wavefunctions compared to those of such single 
band models^. In particular, we found that the low- 
energy quasiparticle band is a result of the crossing be- 
tween the bands of spin-i and spin-| polarons. The | 
polaron has a local spin-1 fluctuation surrounding a spin- 
^ core. The spectral weight for electron removal is ex- 
actly zero at in the region of A: = and (tt, tt) because the 
I band crosses below the | band. The spectral weight 
is also exactly zero at fc = (0, tt) even though there is no 
band crossing. We found that this is due to the orthog- 
onal symmetries of the lowest k = (0, tt) state and the 
state created by removing a k — (0, tt) electron from the 
AFM GS. These findings prompt a more general study 
about more doping holes; here we present results for a 
large cluster with two holes. 

The numerical modeling of multi-hole systems at zero 
temperature is challenging due to the combinatorially 
large Hilbert space and the fermion sign problem. In this 
work, we further validate a recently proposed numerical 
scheme^ as an efficient and systematic way of modeling 
the doped AFM relevant to the lightly doped regime. We 
then use it to calculate the explicit wavefunctions for two 
holes introduced into a half-filled cluster of 32 copper and 
64 oxygen sites with periodic boundary condition. The 
numerical solution points to a competition between local 
attractive potential and kinetic energies. 

The article is organized as follows: in Section II we 
introduce our model starting from a general three-band 
model, and then specify the various assumptions used to 
simplify it. In Section III we provide computational de- 
tails and validate our numerical approach. Section IV 
contains our results. Section V further discusses the re- 
sults and also points out features that are expected to be 
robust in larger systems. Section VI contains the conclu- 
sions. 



II. THE MODEL 

We start from a three-band p—d model which describes 
the basic physics of a hole-doped, charge-transfer gap. 



insulating spin-i antifcrromagneti^ — 

H3B = Tpd + Tpp 4- Aprf 2_^ ni+^^cr 

+Upp ^ ni+^^^ni+^^i + Udd ^ ni,-\ni^i. (1) 

Here, ni,„ = d\,a'^i,'y *^nd ni+^^„ = p\+^_c,Pi+t,c count 
holes with spin a in the Cu id^2_y2 orbital at site /, 
respectively the O 2px/y orbital located at I + ^x/y, and 
Udd > Upp > Apd describe Hubbard and charge-transfer 
interactions. 



Tpd = tpd^[{~pl+,^„ + pl_,^^)di,^ + h.c] 



and 
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PP 2^ ^^Pl+t+S,aPl+<^-<' 



describe Cu-0, respectively 0-0 hopping, where the sign 
ss = ^x^y/\5x^y\ of thc liopping matrix elements is de- 
termined by the relative phases of thc initial and final 
orbitals. The meaning of the vectors e and 5 is explained 
in Fig. 1(a). We have taken the direct Tdd hoping to be 
negligible because of the large spatial separation between 
Cu^+ ions in the lattice. We have also discarded T^^, 
the direct next nearest neighbor 0-0 hopping, which 
was found to have negligible effects on the single polaron 
states found in Rcf. 57] and summarized below in Section 
II. 

At half-filling, in the insulating parent compound, 
states of mainly oxygen 2p character form a completely 
filled band and there is one hole per Cu. If a no-double 
occupancy restriction is enforced, the resulting ground- 
state has AFM order, with a nearest-neighbor Cu-Cu su- 
perexchangc interaction mediated by virtual hopping be- 
tween oxygen 2p and Cu Zd orbitals. An effective model 
for two doping holes, located at O sites, interacting with 
this antiferromagnetic background, can be derived by di- 
rect generalization of the Udd -^ 00 Rayleigh-Schrodinger 
method described in the supplementary material of Ref. 
[571 for the single hole case. If all hopping integrals are set 
to zero in the three-band model (Eq.[T|), thc two-hole GS 
are of the type: 



p\+e,A'- 
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(2) 



with each ket \(Ji")i" specifying the spin of its Cu. These 
states have a degeneracy of (^1^)2^+^. Thc Rayleigh- 
Schrodinger expansion operates in this highly degenerate 
subspacc of the iV-t-2 hole sector of thc Fock space to give 
an effective Hamiltonian for these states. The outcome 
is slightly different than in the single-hole scenario. In 
particular, the projector used in the expansion is 

P2h = J|(l-n^+e,■f'^^+e,;)]^(?^^,t+"■^,;-2n^,^n^,;), (3) 

and allows only states with a full lattice of copper spins 
and no double-occupancies (due to Upp/dd)- In other 



words, the doping holes are forced to hve on O sites. 
After the expansion, the many resulting terms can be 
grouped such that the effective Hamiltonian is written as 
the sum of two parts: 



when close to each other. We write: 
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P-lhHiP-ih 



Ho 



(4) 



Hi is the collection of terms when the two holes do not 
directly affect one another. This single-hole effective 
Hamiltonian was derived in Ref. ISTt 



Hi — Toi 



-t;, 



Hj^i + ffjdi , 



(5) 



where the bare oxygen-oxygen hopping of a hole, Tpp, is 
supplemented by: 

Tswap = -~tsn,^s,^p\^^^,^,^pi+,^c-\a'i^^){ai^J (6) 

^./pd = Jvd 2^ Si ■ Sl±e (7) 

^,/dd = Jdd ^ Sl±2e ■ SlUcril - ni±e^a) (8) 

The physics described by these terms has been discussed 
in some detail in the supplementary materials of Ref. 
ISTI . Briefly, Hj ^ is the exchange between the spin of 
a doping hole located at an O site, with that of its 
neighboring Cu spins. This term favors the formation 
of a three-spin polaron (3SP) )^^'^° either | fl") or | JJ-) ~ 
the corresponding cigcnfunctions for the central Cu-0- 
Cu spins are listed in Table HI These describe the fer- 
romagnetic core of the hole-induced disturbance in the 
otherwise AFM background. Hj^^ is the usual superex- 
change between neighboring Cu spins, which, however, is 
blocked if a hole is located on the ligand O. Such blocking 
decreases the penalty for having the 3SP ferromagnetic 
core in the AFM background. Finally, Tswap describes 
the processes sketched in Fig. [Ijb), where the hole from 
a Cu neighbor to the doping O hole first hops to another 
of its three hole- free O neighbors, followed by the original 
hole falling into the Cu orbital. The opposite phase be- 
tween Tpp and Tswap favors the coherent propagation of 
the three-spin polarons. Even in the single-hole scenario, 
the physics encompassed by these terms leads to qualita- 
tive differences when compared to the t-t'-J modclj^ 

Now we identify the important terms contributing to 
H2, which describes the evolution of the two doping holes 

I 



keeping up to fourth-order terms in the Rayleigh- 
Schrodinger expansion. In the expansion, the appear- 
ance of (1 — P2/1) dictates that all terms must have an 
even power of tpd because the final states must have one 
copper spin per unit cell. Because H2 accounts for only 
2-hole corrections to P2hHiP2h, one can deduce that all 
terms in the second-order 7J, share a factor of tp„. All 

factor be- 



terms in the third-order Hr, share a t? .t 
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cause the second step of any three-step i^^ process would 
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yield virtual excited states projected out by (1 
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TABLE I. Single-hole eigenstates of Hj 



ply creates an oxy- 



gen hole and the arrows in the ket indicate the spins of the 
two copper sites neighboring the oxygen hole. 



The fourth-order H2 collects terms proportional to t. 



PP' 



The second order corrections are Tp- processes that 
link initial and final states in which oxygen holes arc 5 
apart (Fig. la). Because we are considering only the a- 
bonding oxygen orbitals, the virtual excitation is a state 
with double occupancy of an oxygen orbital; that is, the 
inatrix element is non-zero only for singlet correlations: 



H 



(2) 



2t. 



-7f E(-i) 



5{6i_-S2) 



'^a'a ' '-5, 



P'p 



Pl+e+Si,a'Pl+<^+^l,aPl+^+Si+S2,l3'Pl+<^,l3' 



(10) 
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with 



Sa'a ■ S/i'p - T Z^ ^"'^ 



a^P'p-: 



i^x,y,z 



where ct^/q, is an element of the pauli matrix in the i*'' 
direction. 61^2 sums over all generic S values (Fig. [T^). 



d{di ■ 62) = 0[1] if Si and 62 are parallel [perpendicular]. 
The matrix clement is non-zero only for singlct-like con- 
figuration and when the two holes are \6\ ~ -y= apart. 
There are 8 non-zero matrix elements in this situation. 
Two of these correspond to 61 + S2 =0 so the static AFM 



It 

exchange is 2 x -jf^ . There are also 6 other ways for one 
hole to "skip" over the other with such a Heisenberg fac- 
tor. 

There is an abundance of terms in third- and fourth- 
order, which are a sub-set of terms studied in the litera- 
ture of 2-band modclsi^"— To provide a simple physical 
picture, we consider only terms that are greater or equal 

8-16-*'-, 16, 



to Jdd^ which is roughly 
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lb X f^-. f 

81 VV V^ 



2t 



VP 



and A„ 



Unp ^ Qtpp. The basis of this approximation 



is the observation that the dominant short-range physics 
should be already captured by the numerous low-order 
terms with an energy scale of tpp and Tswap,Jpd,Jpp ^ 
0.66f pp, which are already greater than the long-range 
physics at the scale of Jdd ~ 0.2ipp. Adding the rele- 
vant short-range corrections with magnitude greater than 
0.2ipp should then be more than adequate. We also dis- 
cards terms that can be factored into Tpp^ T^wap, Jpd, Jpp 
and the identity processes. These terms effectively scale 
processes by some overall factors which roughly cancel 
out when all parameters are divided by Jdd in order to 
use dimensionless parameters, since Jdd also undergoes 
similar renormalizations in higher order. 

As discussed above, all third-order corrections have 
a prefactor of tppt^^. The perturbation goes through 
two virtual states so the largest possible magnitude is 

I 



^pd ^ pp 



A„rfA„ 



9 PP' 



The splitting due to such a 

pair of Hermitian matrix elements is ~ |ipp, of the same 
order as the Jdd splitting. These processes involve vir- 
tual excitations with no double occupancy. They are all 
multiples of tpp, Tgwap, Jpd and the identity processes and 
arc discarded due to the rescaling argument discussed 
above. Because Upp ~ Apd, other processes would have 
denominators that arc at least a factor of two greater; 
that is, constructive quantum interference is required for 
any terms to be non-negligible compared to the superex- 
change. Constructive interference among tppii^ processes 
requires the same orbital occupancy in the initial and fi- 
nal state, and this can happen only if the two oxygen 
holes are \5\ apart. The transition of interest is then 



p\±e^,aP\±ey,p4,i ^ P'+- -'P^+'= «'^^-" ^^^ effective 



matrix elements have the form {a' , P' , •y' \H2 |a,/3,7), 
a local three-spin ring involving the copper spin sand- 
wiched by the oxygen holes. The correction term can 
be expressed as a summation over each copper spin with 
two additional vector A,j./y summed over ie^/y Not- 
ing that the 3-step hopping would give an overall phase 
of —1 (Fig. [TJi) and that the virtual states have double- 
occupancy, it is not surprising that all operators contain 
a shifted Heisenberg factor: 



H. 



(3) 



^''dp'^PP 



^ pp\^pd 'V ^pp) 



^''dp^PP 

{Apd + Upp)' 



+ [Sa'a-Sl3'p- -^jpl+A,,'yPl+A:,,apl+^^^a,Pl+Ay,p\l,l3'){l,j\ 

+ rSyj-^a'a - -^j pUA^,YPl+A.,'^PUAy,a'Pl+Ay,fs\l, I3){l,j\ 

+ f V/3 -^Yl - ^j p|+A,,,3':P'+A.,a_pJ+A^,yP/+A„,;3|^,a)(^7l 

Y^ Sap ( S'/3'/3 • Sy^ - -j pI+ A^^aPl+A^,ap]+A^^p,pl+ Ay ^p\ln') {I hI 



+ Syl3 iSa'a ■ Sp-p - ^ j p!+A,,T.-Pi+A^,aPl'+A^,Q,'Pi+A^,/3|^,/3')(^7l 



(Sa'a ■ Sy^ - ^jp|+A,,a'Pi+A.,ap|+A„,/3P'+A„,/s|',7'>(^:7l (H) 



r 



These terms are finite if there is a 2-spin singlet 
amongst the 3 spins. The first four terms give an eigen- 



value of +4 



Upp{Ap^ + Upp) 



0.2tpp for oxygen-oxygen sin- 



we can approximate this term by simply raising the en- 
ergy of oxygen-oxygen singlet pairs accordingly. 

In the fourth-order, all tp^ processes can be factored 



glet pair and for triplet pair. The last four terms give 

t^ t 



into two T,, 



an eigenvalue of 
and 0, -3 *'"*'" 



(Api + Upp)-^ 



[Api+UppY^ 



'^swap or Jpd processes. The tz^d^pp processes are 
smaller than Jdd by a factor of ~ 4x4 and no constructive 
interference is possible. 2tp processes are even smaller. 
-QMtpp for triplets. Therefore Therefore, we set h!^^ k 0. In summary, the 2-hole 



0.02fpp for singlet pairs 



correction is in essence: 



H2^ + 4l^Y.^-if^''-''^ (-s^,^-^,,p~\ 



Pl+t+Si,a'Pl-+<^+^i'OtPl+f.+5i+S2,l3'P'-+<^^P 



•-^PP Z^ [^a'a ■ Sp'p ~ 4 ) Pl+e+S,a' 



Pl+c+S.aPi+,^l3'Pl+e,f3 



(12) 



with J, 



(2) 



and J, 



(3) 



PP 



interaction due to second 
Using tpd — 1.3eV, tr 



4f^ 



Upp{Apd+Upp) 

and third order 



i-pp 



0.65el/, A 



pd 



for hole-hole 
corrections. 
3.6eV, and 



U, 



PP 



4eF,— we scale the parameters in units of Jdd to 



find their dimensionless values to be t. 



2.98, Jpd = 2.83, and J, 



(2) 
PP 



PP 



1.3420, and 4p^ 



4.13, tsu, - 
0.9182. 



III. THE COMPUTATION 

Ahhough the thermodynamic limit is important in 
condensed matter theory, models detaihng the interac- 
tions between fermionic carriers and a quantum antifer- 
romagnet present formidable challenges to analytical so- 
lutions. Finite cluster numerical calculations are there- 
fore a valuable tool in extracting information about the 
model of interest. Numerical studies of phase transitions 
would require iV — >■ oo extrapolation; instead, this work 
focuses on the correlations between two holes injected in 
a half-filled A^ = 32 Cu02 cluster with superior topolog- 
ical properties compared to smaller clustersi^ The com- 
putational breakthrough described in this section allows 
us to explicitly obtain the low-energy wavefunctions. We 
can then calculate any correlators of interest, in order to 
understand the nature of these eigenstates. We can also 
perform a direct comparison with the results of the t-t'-J 
model with two holes on the same TV = 32 cluster J^ 



A. An insight about antiferromagnets 

Even after exploiting translational and spin-projection 
symmetries, two holes injected into a half-filled A^ = 32 
Cu02 cluster have a Hilbert space with 0.154 x 10^^ 
states. A system of this size challenges the capability 
of all unbiased methods at zero temperature. One way 
forward is to identify how to drastically reduce this di- 
mension while keeping the most important states in the 
basis. The two doping holes contribute at most a factor 
of 4 ( ^^ ) , so it is the AFM background that is primarily 
responsible for this large number. 

We have previously proposed an efficient numerical 
approach for the modeling of antiferromagnets,— upon 
which we will build here. For completeness, we first 
briefly review the undoped scenario^ - here the system 
is described by AFM Heisenberg superexchange between 



I 

neighbor Cu spins, see Eq. ([5]). If the Cu spins on the 
square lattice are divided into two sublattices, A and B, 
such that spins from each sublattice only couple to those 
of the other sublattice, one of the many measures of the 
staggered magnetization can be expressed in terms of the 
total spin S ^ Sa + Sb- 



- ]^ (E(-i)'^'^'y - ii(2^i 



25'! - S^). 



(13) 

The undoped GS is known to be a singlet, S = 0. For 5*^1 
and Sb to add to zero, their quantum numbers must be 
equal, 5*^ = Sb- Moreover, there are accurate estimates 
of m, which range from ^ 0.3 as A^ — > oo to ~ 0.45 
for A^ = 32ji In order to get such large values for m, 
the sublattice spins S^/b have to be within jg of their 
maximum values. In other words, in each sublattice -^ 
spins add to a total spin of zero while the rest add to the 
maximum 3^. Based on this insight, we showed that 
the GS can be captured systematically by considering a 
subspace specified by a completeness parameter 

Cse [0,1]. (14) 

The subspace equals the full Hilbert space when Cs — 1- 
For smaller values of Cs , only the states with 



S 



A/B 



> 



N 



(1 - Cs) 



(15) 



J is indeed the 



are included. C's = is, therefore, a singlet containing 
the classical Neel state. A linear decrease of Cs yields 
a combinatorial decrease in the number of states. In 
Ref. [5^ we postulated apriorily that the subspace with 
Cs = J captures the essence of the wavefunction, and 
then showed that the convergence is exponential as Cs is 
tuned from zero to unity and that Cs 
"sweet spot" . 

For a large sample, two additional holes cannot drasti- 
cally change the entire AFM background; therefore, the 
stable, systematic convergence of the undoped AFM is 
expected to hold for the two-hole scenario. 

We first test the method against the exactly solvable 
scenario of one hole on a cluster with 32 Cu and 64 
O.— Using the same truncation criterion, we calculate 
the one- hole ground state for increasing Cs- The con- 
vergence is illustrated in Fig. [21 The energy computed 
for Cs = J and ^ is within 3.6% and 0.5%, respectively. 




FIG. 2. Convergence of the one-hole ground state energy. 
Left: GS energy calculated for increasing Cs- The value ap- 
proaches rapidly the exact value marked by the dotted line. 
Right: fractional change in the GS energy for the next incre- 
ment of C3. Solid lines are extrapolation and linear fits. 
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TABLE IL (N=32) The S=| k = (f , f ) single-hole ground 
state's probability in subspaces of particular Sa ® Sb values. 
Numbers are percentages adding up to 100%. Two lines are 
drawn to highlight the Cs = j truncation which discards 
states with Sa,Sb < 6 and yields an energy within 3.6% of 
the exact value. 



of the exact value. Table [TT] shows the exact wavefunc- 
tion's probability in various subspaces Sa ^ Sb of the 
sublattice total spins Sa/b- While the Cs = j trunca- 
tion captures ~ 95% of the probability of the undoped 
wavefunctiou)^ Table |lT] shows that Cs = j captures ^ 
73% of the one-hole groundstate. The next two incre- 
ments of Cs contains 17% and 7%, respectively, of the 
remaining weight. The addition of a hole couples the spin 
background to adjacent values of Sa^Sb, and states 
added by increasing Cs have decreasing importance in 
the wavefunction; therefore, it is reasonable to expect an 
increment of SCs ^ jj from Cs = -j to sufhce. At least 
in the very dilute limit, the truncation scheme provides 
a good starting point to systematically capture the low 
energy state. 

These observations provide reasonable merits for the 
application of this scheme to the two- hole scenario, whose 
ground state should be captured in the subspace of 
Cs ^ j + 0{5Cs) for large N . To provide a conservative 
error analysis, however, we would aim for a capability of 



up to Cs ^ 5, which in turn limits the cluster size to 
A^ = 32. We showed in our previous work^ that a sin- 
gle oxygen hole induces a local disturbance in the AFM 
background and that, for the parameter range of inter- 
est, the disturbance affects ~ 6 — 12 spins around the 
hole. Therefore, a cluster with 32 unit cells should be 
large enough to accommodate two holes without artifi- 
cially forcing them too close together. As discussed be- 
low, our results are different from those yielded by similar 
attempts restricted to A^ = 16 unit cells <^"— 



B. Implementation 

The Hilbert space of the two-hole problem (Eq. |4] to- 
gether with Eqs. Eland [12]) contains (2|^)2^+2 states of 

the form pj+, ,,p|,^,,_„, 11 W'). with / + e 7^ /' + e'. To 
perform computations for large A^, this basis must be 
reorganized to take advantage of translational symme- 
try, total-spin symmetry, total-spin-projection symme- 
try, and, most importantly, the truncation scheme dis- 
cussed in the previous section. The implementation is 
detailed in Appendix A. 

We emphasize that the truncation scheme is applied 
only to the AFM background, i.e. to the parts of the 
wavefunction describing the spins at Cu sites. There is 
no restriction for the two doping holes residing on the O 
sites, apart from no-double occupancy. 



IV. RESULTS 

We apply the approach described above to solve the 
two doping holes problem for the A^ = 32 cluster. We find 
that the low-energy states have a total spin of St = 
or St = 1 (the total spin includes both the contribution 
of the AFM background and of the doping holes). Sec- 
tors with higher St have higher energies, in accordance 
with the trend of finite size AFM computations.— As dis- 
cussed in Ref. [53, the one-hole GS are degenerate at 
k = (±^,±^), so the important total momenta for two 
holes are K = (0,0), (tt, 0), (tt, tt). The convergence of 
the lowest states at these high-symmetry points is shown 
in Fig. [3l The trend of exponential convergence is simi- 
lar to that found for the undoped^ and single hole cases 
(Fig. [2]), signaling that the dominant part of the Hilbert 
space has been captured. We first present the energetics 
then illustrate the numerical eigenvectors in details. 



A. Energetics 

Figure [3] shows the dispersion of the lowest two-hole 
states versus the total momentum K, for St = 0, 1. In 
the St = ^ sector, the K = (tt, tt) global groundstate is 
doubly degenerate. The K = (tt, 0) and K ~ (0, tt) states 
are ~ QAlJdd higher in energy, while the K = (0, 0) state 
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FIG. 3. Convergence of two-hole low energy states at sev- 
eral high-symmetry points. The horizontal dashed and dotted 
lines are the AE^/n = and A_Ejv — levels, respectively, 
i.e. zero binding energy levels calculated at fixed doping con- 
centration and fixed lattice size, respectively. The two values 
bound the binding energy value and converge as Af — !> oo. 



is almost 2 J higher in energy. The shape of the energy 
dispersion mirrors that of some previous studies treating 
iV > 32 systems with two holesJ^i^ Although finite- 
momentum superconducting pair has been proposed by 
Fulde, Ferrell, Larkin, and Ovchinniko v^^i^"^ , this sce- 
nario is not superconducting. For N = 16, the GS is 
instead found to be located at K = (tt, 0).— This signals 
that iV = 16 is definitely too small to capture the two- 
hole physics, which is not surprising since a single polaron 
was found to disturb ^6—12 spins in its vicinityj^ 

As for the St = ^ sector, these states cross below the 
St = states in the region around K = 0. The lowest 
5t = 1,A' = level is doubly degenerate as well. Its 
energy is ~ O.SJdd, equal to the ~ -^ free magnon gap 
of the undoped A^ = 32 clusterji The following section 
confirms that these St = 1, AT = low energy states have 
the two holes in their doubly degenerate global GS with 
St = 0, at = (tt, tt), plus a Q = (tt, tt), S* = 1 "magnon" - 
like excitation in the AFM background. It is thus rea- 
sonable to expect the St = 1,A' = states to become 
degenerate with the K = (tt, tt), 5t =0 states as A^ ^^ oo 
and the free magnon gap closes. Similarly, the low en- 
ergy St = i excitations at K = (tt, 0) and K = (0,7r) 
are the St — two-hole state at K = (0, tt), respectively 
K = (tt, 0), plus a Q = (Tr,-^),^ = 1 "magnon". They 
are also expected to become degenerate as A^ — > oo. 

For the one-hole case, we were able to establish the 
robustness of spin-1 excitations because we found that 
£"3/2 — Ei/2 is much lower than any finite-size spin exci- 



tation, £"„ 



El 



II' 



57 



In the two-hole case, however, 



all spin excitations are higher than Emagnon — Eq, and 
no conclusions can be drawn. 

We note that the two-hole bandwidth is ~ 2 J^d, which 
is roughly equal to the single-polaron bandwidth.— This 
suggests that the GS has very weak, if any, binding. (The 
~ 2Jdd scale has further implication when the wavefunc- 
tion is considered in Sec. IV.) The binding energy is the 
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FIG. 4. Energy of the lowest two-hole states vs total momen- 
tum along high-symmetry cuts, for St = 0,1 and Cs ~ |. 



difference between the two-hole energy and twice the one- 
hole energy, shifted by the energy of the undoped system: 



AE ^ (E; 



2/1 



E, 



Oh) 



2(Ai/i — Eoh). 



(16) 



This value must be extracted with care for a finite cluster. 
Although the A^ = 32 cluster can accommodate two spin 
polarons, and each disturbs ~ 6 — 12 Cu spins in its 
vicinity)^ as we show in the next section we find that 
the low energy states have a low probability to be close 
enough to share a common Cu. Unlike the N -^ 00 limit, 
the local range of this blockade is not negligible in a A^ = 
32 cluster so the kinetic energy is over-estimated. The 
one-hole scenario does not have this artifact, leading to 
an overestimation if the value is evaluated from energies 
calculated at a constant A^, 
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N 
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(32) 
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(32) ^ 



Eoh ) - ^(,^1/1 ~ -^o/i ) 



(17) 



Given that A^-dependent scaling is impossible because 
A^ = 16 clusters have a different ground state (e.g. at 
momentum (0,7r) instead of (7r,7r))^ and that a larger 
cluster cannot currently be solved, we estimate a rea- 
sonable lower bound for the binding energy by consider- 
ing the one- and two-hole energies at fixed concentration 

_n_ _ _2_ _ J_ 
N 32 16' 
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/p(32) 



-,(32) 
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Oh ) 



These definitions converge in the large N limit, 
AE = hm AEn = hm AEil , 

JV-i-oo N^oo " 



(18) 



(19) 



and thus give the upper- and lower-bound for the A^ == 32 
cluster. 

As shown in Fig. |3l the GS energy was extrapolated to 



heE. 



(32) 
2h 



-69.917 Jdd- The AEn = and AEi 







levels in the same figure indicate a binding energy of 
(-0.091 ± 0.272) Jdd- This suggests weak, if any, bind- 
ing, although solutions for larger clusters are required to 
reduce the error bars sufficiently to be able to draw a 
rigorous conclusion. 



B. Wavefunction analysis 

1 . Symmetry 

We define Pxy/xy for reflections about the diagonals. 
The two degenerate S = 0, K = {tt^tt) GSs were found 
to have P^y = ^Pxy = ±1. The two degenerate S = 1, 
K = (0, 0) states related to them by a gapless Q = (tt, tt) 
magnon were also found to have P^y = ~Pxy = ±1- This 
is p-wave parity, but it is not a concern because we are 
outside of the superconducting region. For example, the 
product of two Pxy ~ ^Pxy ~ ±1 wavefunctions would 
yield a dx-2-y2 symmetry. Defining Px/y for reflections 
about the lattice parameter directions, the K = (0,7r) 
state was found to have py symmetry and K = (tt, 0) 
state has px symmetry as in the t-t'-J model. 

Unlike in the A^ = 32 t-t'-J model whose lowest two- 
hole S = 0,K = state has s-symmetrjii^, the lowest 
S = 0,K = state has Pxy = Pxy = —1, which is 
dx^-y-i symmetry. Note that the lowering of s-symmetric 
state in the t-t'-J model is due to the t' hopping of the 
Zhang-Rice singlet, and the t' hopping is included in this 
work due to the explicit consideration of oxygen-oxygen 
hopping. 



2. Charge correlations 

The charge correlation as a function of the separation 
R between the two doping holes can be characterized by 
the expectation value of: 





(20) 


C{R) = Y. c{ij) 


(21) 


\,,\=R 





with '^jiC{R) = 1. On a finite cluster, the number 
of hole-hole configuration separated by a distance R > 
■J is limited by the periodic boundary condition to be 
smaller than in the LF' — N ^f oo limit. Therefore, such 
correlations should be compared to the probability P{R) 
of two randomly distributed oxygen holes to be separated 
by R in the same finite cluster. Then, the correlations 
can be gauged by 



AC(i?) ^ C{R) - P{R), 



(22) 



of a 



meaning that the correlation is the same as that 
random distribution if AC{R) ^ {AC{R)} = 0. 

The expectation values C(R) = {C{R)) are shown in 
Fig.[5]at various high-symmetry points. The lowest states 
with K = (tt, tt) St = and K = (0, 0) S't = 1 have the 
same C{R)^ confirming that they are indeed linked by a 
Q = [it, tt) magnon excitation. 

AC{R) features a monotonic increase with R for the 
GS, a small peak of correlation at _R = 2a for (tt, 0), and 
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FIG. 5. Top: Probability of charge separation, {C{R)) for the 



lowest state in the Cs 



subspace. Bottom: The difference 



from a random distribution. See text for details. 



local attraction for the 5* = A' = (0,0) state. Sec- 
tion IV. A shows that the ~ 2Jdd single-polaron band- 
width is the energy difference between the GS with ex- 
tended AC{R) and the K = (0, 0) dx2_y2 state with local 
AC{R). The low energy two-particle band is thus a re- 
sult of the kinetic energy competing with an attractive 
energy which induces a local dx2_y2 pair. 



3. Polaronic Nature 



In the single-hole case,— at low-energies the mobile 
carrier is well described as the spin-i 3SP's, | ft) and 
I ||.) in Table HI which are correlated states of the oxy- 
gen spin with its two neighboring copper spins. For the 
exact ground state (Hj^^) ^ —0.9Jpd, close to the —Jpd 
energy of the exact 3SP. The two-hole solutions yields 
{Hj ^) ^ —l.SJpd, showing that thinking in terms of 
3SP is still valid and fruitful. In this scenario, the oxy- 
gen holes could neighbor the same copper spin to form a 
5-spin object, but C{R) of Fig. [5] shows that the proba- 
bility for this is low. Ignoring these shared-copper config- 
urations, the wavefunction contains two polarons involv- 
ing a total of six spins. Noting that {Hj ^) ^ —l.SJpd, 
the single-polaron levels in Tab. U suggest that the domi- 
nant part of the wavefunction contains four possible 3SP 

pairs with — p^ 



-2: I ^)| 4), I ^)| 4), I ^)| it), 
and I ^) I 4|) . Other pair configurations would have 

{H., 






£ {0,±^,±l}, and can be ignored. This can be 



achieved using a 3SP-pair projector 



Pssp -= n 

AAG{0,±i,±l} 



Hj^JJpd - AX 
-2-AA ' 



(23) 



wherein terms in the product vanish for the eigenval- 
ues of the excluded levels, while scaling the relevant 
{Hj ^) ~ —2 Jpd levels to unity. P'^sp thus projects out 
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FIG. 6. Top: Probability of singlet 3SP pair, {S3Sp{R)) for 
the lowest states in the Cs ~ \ subspace. Bottom: The 
difference from a randomly distributed paramagnetic config- 
uration. 



the subspace of having two 3SP's with any spin-spin cor- 
relation. All low energy states have > 0.78 probability 
in the P35p-projectcd subspace, roughly the square of 
single-hole solution's ^ 0.9 overlap with a 3SP.— It is 
thus clear that the extra holes form three-spin polarons. 
Section IV. A shows that the low-energy two-hole states 
are total spin S = Q singlets. It is thus interesting 
to check if the wavefunctions can be described simply 
by a S" = AFM background plus a 3SP singlet pair, 
^ (I H")! il-) ^ I tr)l -U-))- In Appendix B, we derive a nu- 
merical operator 



S3Sp{R) = ^ ssspiv), 



(24) 



\ri\=R 



which measures the probability of having such a singlet 
as a function of charge-charge separation R, within the 
subspace projected by P^sp- The value is normalized 
with the analogous triplet measure, defined in a similar 
fashion: 

^(53Sp(i?)+T3Sp(i?)) = l. 

The singlet correlation should be compared to 
J2RP'i'^) ^ Ij th^ random distribution of two 3SPs 
spread over the 64 oxygen sites but with no shared cop- 
per spin. The probability of singlet correlation is -j in a 
paramagnetic state. The difference of interest is then 



A^35p(i?) = S^spiR) 



P'{R) 



(25) 



Figure m shows S^spiR) and AS^spiR) ^^r the high- 
symmetry points eigenstates. It is evident that all states 
have enhanced short-range singlet nature as compared to 
the random distribution - a feature absent from previous 
studies using small clusters or the t-t'-J model which does 
not distinguish oxygen and copper sites. This short-range 
nature is not an artifact of small cluster squeezing the 



two polarons together because the cluster does allow the 
two holes to separate further apart, as confirmed by the 
long-range nature of charge correlation in Fig. [SJ 

ASsspiR) tells only part of the story because it 
is counter intuitive to have the two polarons with in- 
dependent singlet-triplet ratio as indicated by Fig.[6l The 
total spin is S't = so the iV — 4 copper spins not neigh- 
bored by the oxygen charges must carry a spin of 1 when 
the two polarons tend to form a triplet. The 6-spin de- 
scription, which is sufficient to describe the location of 
charge, is not sufficient to capture the spin carried by 
the polarons. The polarons produce a spin disturbance 
with a i?-dependent spatial extent beyond the two cop- 
per spins sandwiching the oxygen charge. Note that the 
non-local spin disturbance around a polaron is already 
observed in the single-hole case.— 



V. DISCUSSION 

Section IV shows that in the GS, the two holes are not 
correlated spatially, since the probability AC{R) to find 
them at i? < 2a is less than for randomly distributed 
holes. This suggests that if the binding energy indeed re- 
mains negative in the limit N ^ oo, the pair would have 
to be correlated in k-space. One needs to be able to study 
larger clusters in order to fully settle this issue. However, 
local pairing correlations are seen for the two holes in the 
S=0, K=0 state, showing that the competition between 
kinetic energy and a local attractive potential. For an 
energy cost of a single polaron's bandwidth, ~ 2 J^j^;, the 
state changes from the K = (tt, tt) delocalized GS to a 
K — dx2-y2 local pair. This strong momentum depen- 
dency suggests that, amidst the presence of a local dx2_y2 
attraction, a local-pairing scenario cannot completely ex- 
plain the full problem, at least within this model in the 
n ~ 2/32 low-doping scenario. 

Although we have broken the technological limit, two 
holes in a 32-Cu02-unit-cell finite lattice cannot rigor- 
ously model all aspects of an infinite system. For ex- 
ample, the possibility of a quantum critical point in the 
4.3% low-doping regime^ cannot be tested. There are 
other questions: what is the implication of the ^ 2Jdd 
penalty when long-range AFM order due to Jdd is de- 
stroyed? Is there a better measure of the binding en- 
ergy than referencing from a single polaron along with 
a AFM background? What are the translational and 
point-group symmetries of the overall wavefunction of 
more than two holes? Finite-size scaling beyond N=:32 
would certainly provide a more robust description of the 
physics but would require more advanced approaches and 
better technologies. 

Next we point out that the weak local attraction, if 
it ever gains prominence as doping is increased, would 
not lead to the unrealistic real-space clustering of holes 
as doping concentration is raised. The only hint of clus- 
tering mechanism is the attraction for local dx2_y2 cor- 
related holes as in the K,St = state. The nature 
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of this attraction can be grasped from the i? < a lo- 
cal peaks of singlet tendency in Fig. [HI By expanding 
the 6-spin wavefunctions 4= (| f|-)| 4) ± | ^)| 1T)) and tak- 
ing the 5' • 5" between a copper spin from the first 3SP 
and one from the second 3SP, the expectation value is 
— joJdd if they are singlet and +jaJdd if they arc triplet. 
Two 3SPs can take advantage of this energy lowering 
by forming a singlet, but only if they are separated by 
an empty oxygen site which mediates an ordinary S ■ S 
Hcisenbcrg AFM bond. The two-hole numerical solu- 
tion indeed shows a local maximum when R/a < 2 
(Fig. [6l). Now, by expanding the 12-spin wavefunctions 
obtained as the direct product of two 3SP singlet pairs: 
1 (I ^)l 4) - I J)l J» (I ^)l 4) - I 4)1 ^)), the expecta- 
tion value of S" • S* between a copper spin from the first 
pair and one from the second pair is exactly zero. There- 
fore, more than two holes cannot lower their energy by 
clustering. 

Lastly, we point out features of the two-hole solution 
that is expected to prevail as the cluster size N is in- 
creased. The doped charge would form a core of 35*^, 
which is surrounded by spin disturbances which deter- 
mine the spin of the polaron, in agreement of the one- hole 
scenario^. The extent of spin disturbance varies with 
hole-hole separation. The R-dcpendent polaron-polaron 
correlations would not change drastically because the 
TV = 32 cluster can accommodate the two holes with- 
out artificially forcing them together as in the A^ = 16 
case. Due to this fact, the observed cross-over between 
locally and non-locally correlated states within a ~ 2Jdd 
energy window is expected to be robust. 



VI. CONCLUSIONS 



To summarize, we have extended our previous 
.Qj.^ ,57,58 ^p study two-hole states in Cu02 planes, in the 



context of the spin-polaron model of Ref. 1571 . Our nu- 
merical approach bypassed various technical limitations 
and extracted the explicit low-energy wavefunctions of 
two holes injected into a half-filled system with 32 cop- 
per and 64 oxygen sites. 

The A^ = 32 solution was found to be different from 
the A^ = 16 solution. Similar to the A = 32 two-hole t- 
t'-J model, the GS was found at AT = (tt, tt), with IS.E ^ 
jT Q = (tt, tt) magnon excitation. The binding energy 
was found to be (—0.091 ± 0.272) Jdd- In contrast to the 
t — t' — J model, we found the lowest state at K = (0, 0), 
without a (5 = (tt, tt) magnon, to have be a dx^^y2 locally 
bound state. 

Further analysis of the wavefunctions revealed that the 
charge carriers are 3SPs, | ■^) and | 4)7 even in this multi- 
hole scenario, but the polarons' spin disturbances are 
extended in range, involving more than the two copper 
spins sandwiching the 3SP. From the correlation values, 
we established that the low-energy band results from the 
competition between kinetic energy and a local attractive 
potential which induces d,j.2_y2 pair. 



Lastly, we showed that real-space hole clustering is un- 
likely at higher concentration and also noted aspects that 
are expected to be robust for larger systems size. Study 
of higher-doping scenario could be interesting but would 
require more advanced technologies. 
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Appendix A: Implementation of the octapartite 
truncation for the AFM background 

First, we define singlet and triplet creation operators 
for the two oxygen holes, when located at sites / 7^ I': 
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Recall that S is the total spin of the A^ Cu spins (the 
AFM background). According to quantum mechanical 
angular momentum addition, upon introducing the two 
doping holes, the A^ -I- 2 spins can add up to a total spin 
St by mixing a two-hole singlet with a background of 
total spin S = St, or by mixing a two- hole triplet with a 
background of total spin S = St, S't ± 1. Taking St = 0, 
any spin background can be specified orthonormally with 
respect to a position I: 



o]Mi 



\a,S' 



0)1 



Ez=-ic(^,^,^T)tw.i'|ai'5'^ 



-^)i 



(A2) 



Here, a denotes a particular group of 25 -I- 1 spin con- 
figurations related by the total spin raising and lowering 
operators 5*^ = ^j Sf^, summed over all Cu sites. The 
total spin of this a group can be 5 = S't for two-hole 
singlet and S = St,St H for two- hole triplets. Due to 
the choice of the overall projection St = 0, the two-hole 
singlet would mix only with backgrounds with S^ = 0. 
The two-hole triplets would mix with three different pro- 
jections \a,S^ = —l)i, \a,S^ = 0)1, and \a,S^ = 1)/ 
from the group a. The weight c{z, S, St) is the Clebsch- 
Gordon coefficients for mixing these three states with the 
three two-hole triplets to achieve a state with total spin 
St and S^ = 0. 

Exploitation of the translational symmetry is per- 
formed by the use of a Fourier series of the form 



E 






'+£' 



I")/; 
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however, care must be taken to ensure orthonormality. 
Due to the commutation relation of the triplet and singlet 
(Eq. IAl|) . the Fourier series is not straightforward for 
e = e' = tr^iy when both oxygen holes occupy x- or y- 
rung oxygen orbitals (see Fig. [1^). The specification of 
these two-hole configurations requires the two orthogonal 
periodic lattice vector Lq and L\ of length a//V for the 
2D, N-unit-cell lattice. Defining 



El = (;[) - lo,l[ - h), 



(A3) 



most hole-hole configurations can be classified in the re- 



gion 



0<dlo< ^,Q<5h < ^ 



(A4) 



These states can be expressed using iV-term Fourier series 



\Oa;x/yy, SI, a, K) = — = ^ ( 



I — h 



.... (^5) 

Because the two oxygen holes are indistinguishable 
fermions (Eq. lAip . there are three remaining 51 val- 
ues which require special attention due to the periodic 
boundary condition. 

5la = ^,5li=Q 



Slo = Q,dli 



hi 

2 



'^^0 = -Y^Sh = Y 



(A6) 



The number of terms in the Fourier series depends on 
the spin background translated by 51: Tsi\a)i. If such 
a translation yields an orthogonal state, i{a\Tsi\a)i — 0, 
the Fourier series still has N terms. For the example of 
5l = [if,0), 



J 



I' — ^xx/yy 1 ^^ 



L -1 M-i 



\a)i 



^°'«^U,„,+.,+..,. (i + -□^■'^"'^r^ ) !«)/, 



(A7) 



r 



where sn is the sign change due to hole-swapping in there can be only -j terms in the Fourier series due to 
the singlet or triplet (D]^,, = soDIJ. For the case of ^he term (l + sae''^°^TLo). For this case the series 
{a\Tsi\a)i = ±1, the above expansion makes clear that >, ■ th f ^ 



^^^-^ 



D,,;yy,5l,a,K) =. J- J2 e^""^'^ E ^'"""'"OU 



l+e^,^,l+Sl+e^,J^)l- 



(A8) 
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The formulation for the case where one hole occupies of 51 arc unique and the Fourier series has the form 
pi+e^: and the other pi+e is straightforward. All values 
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Therefore, a full orthonormal Hilbcrt space 
can be specified by the states \Dxy,dl,a, K) and 
\Dxx/yy,Sl,a, K). Translational symmetry is specified 
by the quantum number K. Total spin and its projection 
are specified by the singlet /triplet nature of the two 
oxygen holes (D) (singlet-triplet) in conjunction to the 
spin background a. 

The above formulation allows the enumeration of dif- 
ferent oxygen-oxygen configurations for the two doping 
holes, for a given |a) background. An enumeration of 
the different \a) states is now the missing step for com- 
putation. This is not a trivial exercise because the com- 
putation for large N requires that a large number of pos- 
sible I a) states be discarded according to the truncation 
criterion of Eq. [151 Note that this criterion is based on 
the total spin of each sublattice, and the basis a is built 
according to quantum mechanical spin addition. 

An arbitrary enumeration of \a) basis would not yield 
an efficient computation because of several challenges. 
First, the truncation of less important states should be 
controllable systematically and flexible enough to adapt 
to the many unknowns of doped systems. Second, within 
the truncated basis transformed away from the natural 
z-projected description, there should be a fast way of 
indexing initial and final states upon Hamiltonian oper- 
ations. Third, one should have apriori knowledge about 
the identically-zero matrix elements in order to avoid 
them for sparse matrices' performance scaling. 

The above formulation specifies a basis state in terms 
of a Fourier sum over a particular reference site I. The 
spin background \a)i and one oxygen hole pj,^ ^ arc ref- 
erenced at I while the other oxygen hole is located at 51 
away from the reference oxygen hole. The hopping terms 
Tpp, Tswap and Jpp operate on both oxygen holes. Hop- 
ping of the non-reference oxygen hole would not change 
the reference point I; however, hopping of the oxygen hole 
at / -|- e would yield a state with no oxygen at / -I- e. This 
simply means that the final state can be specified at a 
new referencing position I' with a phase shift of e*^^('~' ), 
and that the new spin background is obtained by a trans- 
lation \a')i' — Tr{l — l')\oi)i. It is thus advantageous to 
build the |a); spin basis such that the construction is the 
same upon translation. 

The octapartite approach introduced in Ref. |58| sat- 
isfies all of the above requirements. First, the copper 
spins are divided into eight groups as shown in Fig. [T] 
For each group of -^ spins, we start with the z-projected 



representation of 2 s states and transform into a Clebsch- 
Gordan basis |sjv,s^). Each z-projected basis state is 

8 S 

represented using bits of an unsigned integer and each 
Clebsch-Gordan series is stored as a sparse vector. With 
a particular enumeration of |sjv , s^) states, we mix two 

identical enumerations to build \sn_^s\) according to 
Clebsch-Gordan addition. Then with two identical enu- 
merations of |sjv,s^), we build a single enumeration of 

4 4 

|s«,s^). Finally, we can similarly enumerate the over- 




FIG. 7. A'' = 32 cluster divided into eight groups labbeled = 
7. The spins within each group are connected by multiples 
of (2a, ±2a). Spins from a particular group are always in the 
same environment, eg 1 always 4,5,6 and 7 as neighbors. 



all I a) background as |sAr,s^) from the enumeration of 
|s« , s^). This is a recursive enumeration procedure. At 

2 2 

each stage, the states are indexed by a non-negative inte- 
ger, state_index. The state represented by a particular 
index value can be derived from the "parent" basis with 
half the number of spins using, for example, the following 
loop structure: 



J 



sizc_t statc_indcx =0; 

//in general total spin can be odd multiples of 1/2, so I worked with 2S 

const int inin2S=minimum_2_tiines_total_spin ; 

const int max2S=inaximum_2_tinies_total_spin ; 

const int niinSub2S=iTLininiuin_2_times_sub_lattice_spin ; 

const int maxSub2S=niaxiirLUirL_2_timos_sub_latticc_spin ; 



for(int ss=niin2S ; ss<=max2S ; ss-|-=2){//siaf es with higher total spin have higher indices 
for (int ssa=0; ssa<=maxSub2S ; ssa-|-=2){//Zoop over the left ket 
for (int ssb=0; ssb<=maxSub2S ; ssb-|-=2){//ioop over the right ket 
if( Clebsch_Gordan_coefficicnts ==0) continue ; 
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const sizo_t na=number_of_ ( ssa + 1) -blocks_wi t h_t ot al_spin_ssa ; 
const size_t nb=number_of_ ( ssb + 1) -blocks_wit h_t ot al_spin_ssb ; 
for(sizc_t aa=0;aa<na;++aa) { 
for(size_t bb=0;bb<nb;++bb){ 

for(int zz=— ss ; zz<=ss ; zz+=2,++state_index ){//ioop over the ss+1 components of z— projection 
f print f ( stderr , "State %llu has s=%4f and sz=%+4f , product of" , 

state_index ,0.5* ss ,0.5* zz); 
f print f ( stderr , "block %llu/%llu of the s=%4f sector of left ket and " , 

aa , na ,0 . 5 * ssa ) ; 
f print f ( stderr , "block %llu/%llu of the s=%4f sector of right ket" , 
bb , nb ,0 . 5 * ssb ) ; 
}//zz 
}//bb 
}//aa 
}//ssb 
} // ssa 
}//ss 



The benefit of this enumeration is an instantaneous 
"reverse lookup" . When performing Hamiltonian oper- 
ations, one is really interested in the non-zero overlap 
between the outgoing states and those in the orthonor- 
mal basis. The naive way is to compute the dot product 
against all basis states, but this adds a 0{n) layer on top 
everything else and is detrimental in the case of large 
systems. Under the above looping scheme, an increas- 
ing state_index is associated with increasing values of 
ssa, aa, ssb, and bb. Because the index shift for states 
of these combinations is known a priorily, any arbitrary 
ket |(Ta,cr^)|cr(,,cr^), is trivially associated with these four 
indices so states with non-zero overlap arc known im- 
mediately, with an 0(1) reverse lookup operation. Non- 
zero matrix elements can thus be computed efficiently by 
changing ssa, aa, ssb, and bb before the reverse lookup. 
The neighboring pattern in Fig. [7] eases the determina- 
tion of spin backgrounds upon translation. The deter- 
mination of matrix elements is trivially parallelizable. It 
is apparent that truncation according to the criterion of 
Eq.[T5]can be performed at last stage of mixing by simply 
tuning the value of minSub2S. 



Appendix B: Singlet correlator between two 
three-spin polarons 

If we consider only the 6 spins involved in the two 3SPs, 
the projected wavefunction is a superposition of the four 
possible 3SP pairs: 

m) = a j= h 6| fr) I fr) 



+c'- 



V2 

4) + 1^)1 4) 

V2 



d\m^)- 



oxygen spms, p\ ^ and P2 



Defining the oxygen-oxygen singlet operator for the two 

(Bl) 






1,2 



1 

V2 



.t ^t 



t ^t 



Pl,\P2,i Pl,\P2,-\ 



r 



the probability of finding a oxygen-oxygen singlet in j^g) 
is 

{si2^i:2) = la^ + l{b' + c^+<f). (B2) 

Solving for a? and generalizing the expression to dif- 
ferent hole-hole separations r], the singlet nature of the 
3SP pair when at distance r) apart is gauged by 

S3SP{V) = TB v^ -/ no ^ (B3) 



{P3SpE,'Hv')P3Sp) 



where 



Kri)-J2^^ 



^l+e,l+e+ri^l+<^,l+<^+V 



measures the probability of the two O holes to be in a 
singlet at distance rj apart, irrespective of what the Cu 
spins are doing. Thus, s^spiv) is the probability, within 
the projected subspace, of the two oxygen holes and the 
AFM background to cooperate to form a 3SP singlet pair 
separated by a distance rj. The value ranges from zero for 
no singlet nature to unity for pure singlet at this hole- 
hole separation; however, these two extreme values are 
not possible due to, for example, the spatial spreading 
required to lower the kinetic energy of Tpp and Tswap- 
The measure can be summed as a function of hole-hole 
separation 



I';l=-R 



(B4) 



which is normalized with the analogous triplet measure, 
defined in a similar fashion: 

^(^35p(i?)+T3Sp(i?)) = l. 



The singlet correlation should be compared to 
J^pP'i^) — 1' the random distribution of two 3SPs 
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spread over the 64 oxygen sites but with no shared cop- 
per spin. The probability of singlet correlation is 4 in a 



paramagnetic state. The difference of interest is then 

P'{R) 



ASsspiR) = S3Sp{R) - 



(B5) 
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